Complex light: Dynamic phase transitions of a light beam in a nonlinear non-local 

disordered medium 
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The dynamics of several light filaments (spatial optical solitons) propagating in an optically nonlin- 
ear and non-local random medium is investigated using the paradigms of the physics of complexity. 
Cluster formation is interpreted as a dynamic phase transition. A connection with the random 
matrices approach for explaining the vibrational spectra of an ensemble of solitons is pointed out. 
General arguments based on a Brownian dynamics model are validated by the numerical simulation 
of a stochastic partial differential equation system. The results are also relevant for Bose condensed 
gases and plasma physics. 



I. INTRODUCTION 



At low temperature, the dynamics of complex media 
is dominated by the potential energy landscape (PEL), 
i.e. the multi-dimensional surface of the potential energy 
as a function of the molecular coordinates. Q, A 
disordered system sampling different PEL configurations 
undergoes "dynamic phase transitions," perhaps one of 
the most spectacular ideas of the physics of complexity 
0, 0, 0) ■ For the dynamic glassy transition, observed 
in a certain class of (glass-forming) supercooled liquids, 
this "configurational sampling" is at the origin of an in- 
crease of viscosity of several orders of magnitudes, up 
to values comparable to solids (for recent reviews see 
for example 0, [Tol|'). This phenomenon cannot be di- 
rectly ascribed to a purely thermodynamic transition; 
conversely, it is now widely accepted as a dynamic ef- 
fect. During cooling, different PEL regions are visited, 
and the almost abrupt change of viscosity is associated 
to the transition from "saddle-dominated" to "minima- 
dominated" PEL basins. Since phonons are the elemen- 
tary excitations around minima of the PEL, this process 
is also denoted "phonon-saddle" transition, and can be 
realized while keeping fixed the temperature and acting 
on some pararneter, like the particle density or the inter- 
action range. 11] 

It is not difficult to recognize the fundamental charac- 
ter of these ideas, and the fact that they are not limited 
to the specific contexts where they originally developed. 
In this article I will show that a light beam propagat- 
ing in disordered medium (i.e. a medium with negligible 
optical losses whose refractive index is randomly vary- 
ing) may undergo a sort of dynamic phase transition. 
This happens when, due to an optically nonlinear re- 
sponse, multiple filaments are generated. Their number 
and properties depend on the mutual interaction range 
and, in essence, they behaves like molecules in a complex 
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medium, exhibiting dynamic phases. The process can 
described in terms of the appropriately defined "inherent 
structures" and "saddles," in perfect analogy with the 
physics of disordered materials. This qualifies as a sort 
of soft, or complex, light. 

Here the word "filament" (roughly) identifies a spa- 
tial soliton (SS), which is a non-diffracting light beam 
generated in an optically nonlinear mediurn, with an in- 
tensity dependent refractive index. [12, [ijj Ujl If a suffi- 
ciently intense laser light propagates in such a material, 
self-induced trapping counteracts the natural tendency 
to diffract, and a tightly focused SS can be observed. 
For example, in nematic liquid crystals (NLC), with laser 
wavelength in the near infrared, it is possible to gener- 
ate very thin (few microns waist) SSs, able to propagate 
undistorted for millimeters: hundreds of times the dis- 
tances attainable in absence of a nonlinear self-action. 

El El III 

At each SS is associated an optically induced perturba- 
tion An to the refractive index, which is at the origin of 
the self-focusing. If the incident beam is sufficiently wide 
many SSs are generated by the same input, as shown 
for example in |l6l ITsl IT^ . These filaments may prop- 
agate with various degrees of interaction, relying on the 
material properties and, in particular, on the so-called 
"non-locality" , which can be kept in mind as the ratio 
between the spatial extension of the refractive index per- 
turbation induced by one filament and its transverse in- 
tensity waist. In NLC, the degree of nonlocality can be 
simply controlled by a voltage bias,f20j and hence the dy- 
namic phase transition may be induced accordingly. ^^Ij 

The following results not only are well suited to de- 
scribe light-soft-matter interaction (as in 0, 0, |23, l2ll 
12^'). but can be applied to interpret the dynamics of ul- 
trashort laser pulses propagating in air, |3 as well 
solitons in photorefractives. 1251 1261 media with thermal 
nonlinearities and plasmas, [271 [23, [^^semiconductors 
[sof , discrete solitons (see for examplepTj and references 
therein), plasmas [2^ and Bose-condensed gases [33.[33j . 

This article is aimed to the introduction of the leading 
idea, and to a phenomenological description, supported 
by Brownian dynamics and stochastic partial differen- 
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tial equations (PDEs) numerical simulations. In order to 
simplify as much as possible the presentation no attempt 
will be made to a theoretical analysis, which is deferred to 
future publications. In section II the model for the prop- 
agation of several light filaments in a random medium is 
linked with Brownian dynamics. In section III evidences 
of phase transitions from numerical simulations are re- 
ported. Section IV is dedicated to the definition and 
the analysis of the so-called inherent structures. In sec- 
tion V the generalized inherent structures are considered 
for the final settlement of the dynamic phase transition. 
In section VI the noise quenching process is addressed. 
In section VII a connection with the theory of random 
matrices and the so-called dynamic structure factor is 
established. In section VIII numerical simulations of a 
well known model in nonlocal soliton theory validate the 
general arguments of the manuscript. Conclusions are 
drawn in section IX. 



II. REDUCTION TO A BROWNIAN 
DYNAMICS MODEL 

For the sake of simplicity, the analysis is done with 
reference to one-dimensional (1+lD propagation) beams. 
Indeed, differently from standard thermodynamic transi- 
tions, dynamic transitions can be obtained in low dimen- 
sional systems. Additionally, the considered case reflects 
typical experimental geometries for the investigation of 
modulational instability (see e.g. 0,^]). The gener- 
alization to higher dimensional problems of what follow 
can be readily imagined. Consider the Fock-Leontovich 
equation, which describes the paraxial optical propaga- 
tion in a nonlinear medium: 



2ik 



OA d^A 



dz dx^ 
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A = 0. 



(1) 



A{x,z) is the complex amplitude of the optical field such 
that / — \A{x,z)\'^ is the intensity, k — 27rn/A is the 
wave-vector and n the refractive index at wavelength A 
in absence of non-linear effects. An is the optical in- 
duced perturbation to n, its functional relation with / is 
non-local: the specific distribution of An(x, z) depends 
on the whole profile I{x,z). Various models are given 
in the literatures relating An and /, as the Kukhtarev 
equations for photorefractives, (for a review, see Del Re 
and coworkers in |l3ll) the heat equation for thermal non- 
linearities. |28l l29j | re-orientational equations for li quid 
crystals, [sj mode coup ling the ory for soft-matter, |23| 
or generic nonlocalities. [S^, |3a, |33, 

For relatively low intensities An is linear with / and, 
in the presence of many incoherent filaments, it can be 
written as the sum of the intensity distributions of the 
filaments, with p = 1,2, ...,N. The overall An is thus the 
sum of the index perturbations of each SS: 



An] I] 



N 

^An[/,] 

p=i 



(2) 



Eq. (0) is valid whenever An/n << 1, as typically veri- 
fied in the reported experiments, and when the filaments 
are mutually incoherent, i.e. the relative phase between 
pairs of them is randomly varying. The latter hypothe- 
sis typically holds in soft-matter, by taking into account 
that the SSs spontaneously generate from noisy intensity 
perturbations, and propagate in a thermally fluctuating 
medium. Furthermore, in the presence of many filaments 
the effects of the relative phases and, more in general, of 
the specific form of An[I] can be negligible. This is anal- 
ogous to the typical approach in statistical physics where, 
quite often, the particular profile of pair-interaction po- 
tentials can be replaced by s ome simple model like the 
Lennard-Jones potential. |33 

Following a perturbative analysis, the effect of the in- 
dex perturbation due to all SSs on the trajectory of 
the generic filament p is considered. For N identical 
stable solitons with intensity bell-shaped profile I six), 
Ans{x) — An[Is]{x), and average position Xp{z), using 
the Ehrenfest's theorem of standard quantum mechanics, 
applied to the Schrodinger-like Eq. one has for the 
generic filament 

dAn/n 

Is(x-Xp) — — — dx, (3) 
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with m = J Is{x)dx the power (per unit length along 
the y direction) into each filament, which plays the role 
of the particle mass. The index perturbation is 
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which used in Q yields 
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with 



Vix) = -- r Ans{i + ^)Isi^ ~ ^)dt (6) 
n J-oo 2 2 

the pair-interaction potential. Many derivations of simi- 
lar results can be found in the literature on solitons and 
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solitary waves (see e.g. |40|, |4l| and references therein). 
Here the analysis is specialized for the potential energy 
landscape interpretation. 

In lO the self-interaction term (p = q) can be retained 
since it clearly gives a vanishing contribution. Finally 



9$ 



(7) 



and $ — ^{xi,X2, ...jXn) is the overall potential energy 
surface (the PEL) given by the sum of pair- wise interac- 
tion terms: 



N N 



■ Xk). 



(8) 



]=i k=i 



The dynamics along the direction of propagation is hence 
formally reduced to an ensemble of particles, evolving 
with "time" z. The fluctuations of the medium result into 
a random contribution to An that can be phenomenolog- 
ically included in the model as a Langevin force rip(z): 



' dz^ 



9$ 

dXn 



-T]p{z). 



(9) 



In the following, I will take for rjp a normally distributed 
white noise: 



< Vp{z)Vqiz') >= SpSpqS{z - Z') 



(10) 



with Sp the noise "power" and the brackets denoting a 
statistical average over disorder. ?7p(z) takes into account 
the fluctuations of the refractive index and defines a re- 
alization of the random soft-medium. 

Note that, in typical Langevin models, the random 
term is accompanied by a dissipative term, which is 
in general dependent on the lossy mechanisms in the 
medium, like viscosity. In the limit of small losses and 
small noise such a term can be neglected, as it will be 
done in the following in order to leave the treatment as 
general and simple as possible. The qualitative agree- 
ment with numerical results in section VIII, and experi- 
ments pTj . supports this approach. 

The explicit shapes of Is and Ang are due to the par- 
ticular nonlinear mechanism; for the present purpose a 
Gaussian ansatz for both of them is appropriate, since 
non-local optically nonlinear media are being considered. 
[Mill 113 Taking 



X 

2^ 



gives 



Is{x) = /oexp 



Vix) = Vo 



(11) 



Ans{x) = Ano exp ( 



(12) 



with Vo = Ano/o[27ri;2ii;2/(i;2_|_^2-)ji/2 ^^^^ ^2 ^ y'^ + w'^^ 
V{x) is a Gaussian function within an arbitrary additive 
constant. It is written as in (|12|l in order to have a van- 
ishing $ when all the solitons are in the same position 
("condensed phase"), is the sum of the variances of Is 
and Ans and provides a measure of the interaction range 
between the SSs for a fixed w, like u/w that will be used 
in the following. Each filament increases the refractive 
index (the medium is assumed to be focusing, Ano > 0) 
and the interaction is purely attractive, so that Vq > 0. 

The previous formulation points the connections with 
statistical physics while reducing the model to a system of 
interacting classical particles undergoing Brownian mo- 
tion, a typical model for colloids, where the formation 
of clusters at the glass-transition is a well-known process 
(see e.g. m^M.M)- 

However, before proceeding, it is fruitful to point out 
a subtle issue associated to Eqs. © and H12I) . It is well 
expected that a finite number of classical particles, in- 
teracting by a purely attractive potential, will oscillate 
around the center of mass at equilibrium (i.e. for long 
times). This implies that only one energy minimum does 
exist and corresponds to the condensed phase. At a first 
glance, no dynamic phase transition, due to local min- 
ima of the PEL, is expected. Nevertheless, in the case 
under consideration, it is clearly not possible to consider 
an arbitrarily large propagation distance (corresponding 
to long times). If losses are negligible and if the max- 
imum observation range (in z) is limited by the spatial 
extension of the sample, damping mechanisms can be ne- 
glected and local PEL minima play a role. 



III. NUMERICAL SIMULATIONS OF THE 
BROWNIAN DYNAMICS MODEL 

A typical distance between resolvable (by means of the 
observation of scattered light from the sample) spatial 
solitons, generated for example by modulational insta- 
bility, can be taken as 6 times their waist w. Hence, 
considering thin solitons with waist w — 5fj,m and tak- 
ing tens of filaments implies an input waist of the order 
of 500/im, which is comparable to those typically em- 
ployed in experiments (see e.g. 0,0|). In the following 
two representative cases will be considered: = 10 and 
A^ = 30 (the principal dynamic phase transition, con- 
sidered below, is obtained up to the largest considered 
A^ = 100, not reported). In all the simulations the fil- 
aments are chosen uniformly distributed at z = 0, with 
mutual distance 6w. 

The noise "power" is measured by the corresponding 
adimensional quantity 



W q2 
^.3/2 1/2 P' 



(13) 



which is taken independent from p, for the sake of sim- 



plicity: 



The amount of noise to be included 



in the simulations clearly relies on the specific material 
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FIG. 1: Filaments trajectories Vs the normalized propagation 
coordinate t for a given noise realization and various values 
of the interaction range u/w (here A'^ — 10). 



FIG. 2: Filaments trajectories Vs the normalized propagation 
coordinate t for a given noise realization and various values 
of the int 



(in particular, on the thermal coefhcient of the refrac- 
tive index and on the sample temperature); however it 
is found that the numerical results are very robust with 
respect to noise, and very similar findings are obtained 
when v is varied by order of magnitudes. Hence only the 
case f — 0.001 will be reported as a representative ex- 
ample. The stochastic ordinary differential equations © 
are solved by a second order scheme, whose accuracy has 
been thoroughly investigated and compared with other 
approaches. |43 The results have been validated by halv- 
ing the integration step and doubling the number of re- 
alization in many cases. 

As discussed above, a dynamic transition is attained 
while increasing the density, or equivalently the interac- 
tion length u/w. 11] Here, this corresponds to increase 
the degree of nonlocality. In figure ^ some realizations 
of the SSs trajectories, obtained by the numerical solu- 
tion of Eqs. © when N — 10, are shown for various 
u/w. The adimcnsional "time" t — z/[w{Va/my^^] is 
used on the horizontal axis. For small u/w the filaments 
propagate in the presence of a reduced interaction. Con- 
versely, while increasing u/w, various clusters are formed 
and their number and positions vary with each realiza- 
tion of the noise. A similar result is obtained in the case 
iV 30 (fig. EJ 

In figures 13 and 0] the "final" (i.e. at a fixed t = 
tmax) position of each SS is shown Vs u/w, with the 
results for 10 noise realizations superimposed. Clearly, in 
certain ranges of the control parameter u/w the statistics 
of the final positions are highly peaked around 2 or 3 
clusters, while they spread over a broad region in other 
ranges. The appearance of an interval for u/w where 
two dominant clusters are generated is evident. This is 
referred to as the "principal dynamics phase transition" . 
In the case = 10 a phase with 3 clusters is also present 
and it is somehow more noisy in the case for = 30. For 
a very large u/w the nonlocality is such that all the SSs 
oscillate around an equilibrium position, this has been 
above indicated as the condensed phase. 

Similar results are obtained for an odd number of fila- 
ments (e.g. N — 21). In that case, in correspondence of 
the principal dynamic phase transition, an additional SS 
is found at the middle of the two clusters. 

An open issue is the existence of fractal structures as 
those investigated in |4^ . 
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FIG. 3: (Color online) Crosses, filament positions at tmax = 
1000 for 10 noise realizations; thick black line, average posi- 
tion (see text); white line, inherent structures (here A'^ = 10). 



IV. THE INHERENT STRUCTURE 

Figures 13 and 0] also show the average filament posi- 
tions Vs u/w, when N = 10 and A^ = 30 respectively. 
Note that the positions are determined at a very large t 
so that the clusters are "stabilized," and that the average 
positions for each filament (thick black line in the figures) 
are shown superimposed (i.e. they are not the average 
positions among all the filaments) , so that the thick black 
line provides a visualization of cluster distribution. 

Disordered phases are alternated with others in which 
a fixed number of clusters is obtained. In order to ad- 
dress the existence of some kind of phase transition it 
is necessary to introduce a "control parameter" , as out- 
lined in the mentioned literature. 0, ^3 With this aim, 
I start pointing out the inherent structure (IS) associ- 
ated to the numerical simulations. Once fixed a max- 
imum value for the time tmax, the final distribution of 
filaments is used as guess for a conjugate gradient min- 
imization procedure that finds the nearest minimum of 
the interaction potential $. The corresponding vector of 
positions {xi, X2, xn) is the IS. 0,1131 Its role is ev- 
ident when superimposing the plots of the average final 
positions for the considered realizations (thick black line 
in Figs. 13 and 21) and those of the average IS (thin white 
line in Figs. 13 and 01). Clearly, the latter provide infor- 
mation on the number and the positions of the generated 
clusters. 

According to the literature about complex media, 
8, 10] the average potential energy e/5 of the IS is an 
appropriate control parameter for the dynamic phase- 
transition. In figures 13 and 13 e/s (in units of Vb) is plot- 
ted Vs u/w. The minimum for the potential energy is ob- 
tained at large u/w when all the solitons are in the same 
position and corresponds to $ = (condensed phase); 
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FIG. 4: (Color online) Crosses, filament positions at tmax = 
6000 for 10 noise realization; thick black line, average position 
(see text); white line, inherent structures (here N = 30). 
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FIG. 5: Average potential energy <E> of the inherent structure 
eis in units of Vo Vs u/w; 1000 noise realizations have been 
considered {N = 10). 



conversely when the solitons are uniformly distributed 
(small u/w) $ is at maximum. Hence, while increasing 
u/w, $ is reduced, due to the coagulation mechanism. 
The data in figures and show a decrease of ejs Vs 
u/w up to the first plateau, at the formation of two clus- 
ters. The corresponding value of ejs/Vo is obtained af- 
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FIG. 6: Average potential energy $ of the inherent structure 
e/s in units of Vo Vs u/w; 100 noise realizations have been 
considered (A'' = 30). 



ter by observing that, when two clusters of N/2 SSs 
are formed (for N even), ^/Vq V^/4 as their mutual 
distance goes to infinity. The plateau for V = 10 corre- 
sponds to e/s/Vb = 25, and to e/g/Vo — 225 for N — 30. 

The trend is conserved for various numbers of fila- 
ments. The scale on the u/w axis changes with N (see 
Figs. Elandini) because the degree of nonlocality needed 
for the phase transition obviously increases with the num- 
ber of filaments. 



V. THE GENERALIZED INHERENT 
STRUCTURE 

For the definitive settlement of the phonon-saddle tran- 
sition, the so-called generalized inherent structure (GIS) 
must also be taken into account. 6, 7, 48, 49] It is defined 
as the nearest stationary point of the PEL (where all the 
forces are zero) to the final configuration. The latter is 
used as a guess in a nonlinear solver (I used the c05pbf 
NAG routine, Mark 19) for the N equations d^/dxp ~ 0, 
whose solution, given by a vector {xi,X2, ...^Xn) is just 
the GIS. The saddle-order Kqis of the GIS is the num- 
ber of the negative eigenvalues (imaginary frequencies) 
of the Hessian of $, [see Eq. H14|l below] calculated at 
the GIS; if Kqis = the GIS is a minimum. 

For any realization of the system there is one IS and 
one GIS. In a phonon-dominated phase the two structures 
are the same and Kqis = 0. Conversely, in a saddle- 
dominated phase Kqis > 0, but it tends approximately 
to zero (on average over many realizations) in correspon- 
dence of the dynamic phase-transition, \^ US 

To understand the physical meaning of the GIS, con- 
sider the principal dynamic phase transition. While the 
IS corresponds to the filaments equally distributed be- 
tween the two clusters, the GIS differs for some of the SSs 
positioned at intermediate places; it denotes the way the 
system may escape from the energy minimum. Clearly, 
if the noise-averaged Kqis is high the probability to find 
a direction in PEL to get out from the local minimum 
of $ is high: it somehow measures the number of es- 
cape directions from the PEL minimum. Actually, the 
average Kqis never reaches the zero as discussed in |49| . 
because this would correspond to a complete freezing of 
the system. 

Figure [7| shows the noise-averaged Kqis Vs u/w when 
= 10. It clearly reveals a saddle-phonon transition 
in proximity of u/w = 8; consistently with the phase- 
diagram in figure |S1 For larger values of u/w, the ad- 
ditional dynamic phase transitions are not well defined 
(i.e. Kqis stays around one or two units) due to the 
limited number of particles. This also clarifies the reason 
for introducing a "principal" dynamic phase transition, 
as done above. The latter, in the considered numerical 
simulations, always corresponds to the formation of two 
clusters. This is confirmed in the case TV = 30 (shown in 
figurelHl) where, due to an increased number of degrees of 
freedom, the transition is more evident, and happens, as 
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FIG. 7: Average saddle-order for the case A'^ = 10, other 
parameters as in fig. |3 
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FIG. 8: Average saddle-order for the case N = 30, other 
parameters as in fig. |S| 



before, when two clusters are generated. The transition 
becomes more evident as larger values of N are consid- 
ered, and it is found up to the largest considered value 
(i.e. N = 100). Nevertheless, the case = 10 shows evi- 
dence of this phenomenon, which is hence observable even 
with a limited number of filaments. Roughly speaking a 
given ensemble of filaments can self-organize in various 
way in order to form clusters. Two symmetrical clusters 
is obviously a strong "attractor" for the system, because 
there is only one way to organize it (conversely a larger 
number of clusters can be formed with different aggrega- 
tions of solitons) . This mechanism strongly resembles the 
cluster formation in coupled chaotic maps, a well known 
example of complex system. [50l | 



VI. NOISE QUENCHING 

The standard deviation Ae/s of ejs has a universal 
trend (with respect to number of filaments, their initial 
distance and the amount of noise) , which is shown in fig- 
ures El and El Moving towards a non-local region (i.e. 
increasing the interaction range u/w) Ae/s grows (in 
the landscape dominated phase [l^l) up to the dynamic 
phase-transition, where small values are again achieved 
(in the cluster phase). This quantity has the same trend 
of the corresponding one investigated, for example, in 
[Slf in a Lennard-Jones material glass. Very similar re- 
sults are obtained when simulating a larger number of 
filaments, and for various values of the noise power. 

A reduction of the noise ("quenching") in correspon- 
dence of the glassy-phase is hence evident. To confirm 
this effect, I show in figures [TTI and [T^ the relative max- 
imum deviation from the average position. This quan- 
tity, denoted Ex is calculated by taking the maximum 
deviation from the average (over the considered noise re- 
alizations) position (x) for each filament (at t = tmax), 
dividing by (x) and then averaging the resulting quantity 
over all the A^ filaments. It measures the noise in the SS 
positions, and reproduces the same trend of Ae/5. Be- 
fore the phase transition all the SSs diffuse into a wide 
region, while after the phase transition they are locked 
inside each cluster. 

In the experiments this phenomenon is resolved in 
time, while changing the control parameter for the non- 
locality (e.g. the voltage in NLC experiments [13 )■ This 
means that, since the medium is fluctuating, before the 
transition the number and the positions of the clusters 
are rapidly varying. Conversely, when the two clusters 
are formed, the intensity profile "slows down" , and the 
noise is quenched; this resembles the "critical slowing 
down" 9] in glassy material system (see also the section 
VIII and the movie cited below). 



VII. VIBRATIONAL SPECTRA AND THE 
RANDOM MATRICES 

A variety of issues spontaneously rises, once some anal- 
ogy with a disordered medium has been ascertained. In 
particular, those concerning the spectrum of the fluctu- 
ations and ultimately the propagation of "sound-waves" 
(or better "displacement- waves" ) . These are associated 
to the vibrational spectra of the ensemble of solitons. 

The analysis of sound-waves is one of the most im- 
portant issues in the physics of glassy systems (see for 
example [s^l and references therein). According to some 
authors, the appearance of ultra-high frequency sound 
can be related to the vibrational spectrum of the material 
and, in particular, and to an excess of states denoted "bo- 
son peak" (see |Tl| and references therein). The trans- 
position of these ideas to nonlinear optical propagation 
is beyond the scope of this article. However, it interest- 
ing to observe that one of the most successful theories of 
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FIG. 9; Standard deviation of the energy of the inherent 
structures for 100 noise realizations Vs u/w (N — 10, t,nax ~ 
1000). 
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FIG. 10: Standard deviation of the energy of the inherent 
structures for 100 noise realizations Vs u/w (N — 30, t,nax ~ 
6000). 

structural glasses, the random- matrices approach, |3.l5^ 
also seems well suited to describe the vibrational spectra 
of the positions of a number of optical spatial solitons in 
a disordered medium. 

Given some configuration of the filaments 
(xi, a;2, xat), which can be either the instanta- 
neous distribution at t„iax, or the IS, [ll| the vibrational 
spectra can be found as the eigenvalues of the Hessian 
matrix Hpq {p,q ~ 1...N): 

^ "5 — 5 — ^ ^pq/ ^ V"{xp - Xk) - V"{xp - Xq), 
ox„ox„ ^-^ 

^ ^ k—\ 

(14) 

with V" = (PV/dx^. Due to noise (and eventually to 
chaos), the considered configuration has a certain statis- 
tical distribution; the problem is hence reduced to find 
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FIG. 11; Maximum relative deviation from the average posi- 
tion, averaged over all the = 10 filaments Vs u/w. Param- 
eters as in fig. |51 
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FIG. 12; Maximum relative deviation from the average posi- 
tion, averaged over all the = 30 filaments Vs u/w. Param- 
eters as in fig. ^1 

the corresponding statistical distribution of the eigenval- 
ues of Hpq . Various approaches have been developed and 
successfull y ap plied to explain some material glass fea- 
tures. l53l | What follows suggests that the random 
matrices approach could be very fruitful even in this field 
of research. 

Consider an experiment in which the position of each 
SS can be retrieved by the scattered light from the top 
of the samp le ( as those in the mentioned literature on 
nematicons |l5l Il6j|): the resulting images appear as a 
number of superimposed fight filaments, each with its 
own few-microns waist, with overall intensity distribu- 
tion I{x,z). If a soliton profile Is is associated to each 
SS, as discussed above, /(x, z) can be interpreted as a 
coarse-grained density of particles evolving along z. The 
squared modulus of the double-Fourier transform of the 
image I{x,z), denoted S{kx,kz) (averaged over a given 
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FIG. 13: Pseudo-color plot of the squared modulus of the 
Fourier transform of the course grained density (overall inten- 
sity profiIe)when A'' = 10, u/w — 1 and tmax = 200, averaged 
over 100 noise realizations. 

number of noise realizations) can be interpreted as the so- 
called "dynamic structure factor" js^ of the soft-medium 
realized by the SSs, which play the role of interacting 
Brownian molecules. S{kx,kz) gives the frequency con- 
tent in kz of the z-evolution of the spatial "mode" at k^ ■ 
In other words, once fixed kx, S{kx,kz) provides infor- 
mation on the dynamics of each intensity perturbation 
with period 2-n jkx- [5^ 

Consider, for example, the case in which all the soli- 
tons travel approximately parallel, with a reduced inter- 
action. In this case, if Ax is the average mutual distance, 
S{kx, kz) is expected to be approximately given by a se- 
ries of peaks around kx = ^nmj IS.x and kz = 0, with 
TO = 0,1,2,.... That is the fcz -bandwidth of each "mode" 
at kx is very small. This happens, for example in the 
case N = 10, when u/w — 1, as shown in the inset in 
figure n Taking the corresponding numerical solution of 
eqs. ©, associating to each trajectory a soliton profile 
(a Gaussian profile in units such that w = \), evaluating 
the squared modulus of the Fourier transform of the re- 
sulting I{x, z), and finally averaging over a given number 
of noise realizations, Fig. E|is obtained, which appears 
as anticipated. 

The aim of the random matrices approach is to predict 
the shape of S{kx, kz) with respect to some control pa- 
rameter. In particular, under very general hypotheses, it 
has been shown that, when the interaction range grows 
(or equivalently the density of the particles is increased) , 
S{kx,kz) develops a Brillouin peak with position that 
is linearly dependent on A:a;.|53| This corresponds to an 
X-shape in the two dimensional level-plot of S{kx,kz). 
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FIG. 14: Pseudo-color plot of the squared modulus of the 
Fourier transform of the course grained density (overall inten- 
sity profile)when A'^ = 10, u/w — 8 and tmax ~ 200, averaged 
over 100 noise realizations. 

Repeating the analysis of figure IT!^ for u/w at the dy- 
namic phase transition, should provide some evidence of 
the predicted X-shape. This is exactly what happens, as 
shown in figure ICT 

The formation of a Brillouin peak is obviously related 
to the geometry of the filaments distribution. However, 
it can be also associated to a movement of energy (or 
equivalently of "mass" , which corresponds to the beam 
power into each filament, as discussed above). Indeed, 
when two clusters are formed and move far apart each 
other, there is an evident energy transfer along x. In 
some sense the cluster movement is related to some wave- 
packet motion, in perfect analogy with sound waves in 
disordered media. 



VIII. NUMERICAL SIMULATIONS OF A 
STOCHASTIC PDE MODEL 

In this section I consider a specific example in order to 
validate the previous arguments. Some numerical simu- 
lations of a stochastic PDE system (which is well known 
in the deterministic limit and has been successfully com- 
pared with experiments) are reported and reproduce the 
features previously described by the Brownian dynamics 
approach. For the sake of clearness and compactness, I 
will show sample results for nonlocalities up to the first 
dynamic transition (i.e. the formation of two dominant 
clusters). Additional data will be reported elsewhere. 

The model equations are those of the so-called 1-l-lD 
exponential non-locality, which well describes solitons in 
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FIG. 15: Three realizations of the numerical solution of the 
stochastic PDE JTHJ for cr^ = 0.12. 



FIG. 17: Three realizations of the numerical solution of the 
stochastic PDE (EHJ for cr^ = 1. 
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FIG. 16: Three realizations of the numerical solution of the 
stochastic PDE lO for cr^ = 0.45. 
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as well as thermal nonlinearities 
|29j, and has been thoroughly studied 
(in absence of noise) in the literature (see e.g. [s^, HI] 
and references therein). Using adimensional variables the 
PDE system (corresponding to Eq. J^l) reads as 



pip = 

--\i;\^ + A7j{^,C)- 



(15) 



In H15|l ^ is the normalized transverse coordinate, C is 
the normalized propagation distance, ip is a complex 
field, corresponding to the electromagnetic field; p is 
the medium disturbance, which can be the director an- 
gle of NLC, or the temperature, or the density of the 
medium, depending on the specific physical system; its 
fluctuations are taken into account by a Langevin term, 
which is a white Gaussian stochastic process, such that 
ivi^, OviC, O) = Hi - OHC - C), and A measures the 
amount of noise, cr^ is takes into account the nonlocal- 
ity; as (T^ = 0, the model reduces to the local integrable 
nonlinear Schrodinger equation; the degree of nonlocal- 
ity increases with as cr^. Equations IjlSfl are solved by a 
pseudospectral approach (see e.g. i54i), which maps ((T5|l . 
via discrete Fourier transform, into coupled stochastic 
ordinary e quat ions, which are then solved by the Heun 
algorithm. [53 

In figures I15I16I and El I show some realizations ob- 
tained with different values of the nonlocality parameters 

. The input profile is taken to be a super- Gaussian: 



V'(e,0)-exp(-(e/10)4), 



(16) 



such that it corresponds to a flat intensity proflle, on 
which various solitons are generated via modulational in- 



FIG. 18: Average intensity distribution at = 4 over 100 
noise realizations as a function of ^ and a^. 




FIG. 19: Maximum relative deviation from the average profile 
versus the nonlocality parameter calculated at ^ = 4. 
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FIG. 20: Two dimensional Fourier transform of the intensity 
profile, average over 100 realization, when — 1. 
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stability and interact while propagating along the C di- 
rection. (The reported results correspond to A = 10"'^). 

Consistently with the analysis reported above, for 
small nonlocality various solitons are generated and 
travel almost independently, while at higher nonlocal- 
ity clusters are formed, and two dominant aggregates of 
solitons are clearly evident when = I. 

In figure El I show the intensity profile a,t C, = A aver- 
aged over 100 realizations, it clearly reproduces the fea- 
tures in figures 13 and 01 where two dominant clusters 
appear at a threshold value of nonlocality. 

Note from figs. I15I17I that, for small cr^, the various 
shots are very different from each other, while at high 
nonlocality the two clusters distribution is "freezed." |^ 

The noise quenching mechanism is quantified in figure 
[T^ where I show the noise figure s = e(C = 4), deter- 
mined as 

and corresponding to the maximum (along the ^ direc- 
tion) relative intensity deviation with respect to the aver- 
age profile at C = 4. Figure [T^ shows a drastic reduction 
of noise at high nonlocality, in correspondence of the two 
clusters. 

Finally, the vibrational spectrum (as described in the 
previous section) is calculated: the intensity distribution 
V'(CiC)P: for a fixed cr^, is Fourier transformed and av- 
eraged over the considered 100 realizations. The result 
for CT^ = 1 is shown in figure EUl and clearly displays the 
X-shapc, addressed above. A similar picture is obtained 
for other values of nonlocality in correspondence of the 
two clusters formation. 



IX. CONCLUSION 

The basic aim of this manuscript is to point out one of 
the possible connections between the physics of complex 



media and that of intense laser light interacting with mat- 
ter. During nonlinear optical phenomena in disordered 
media (or even in media where the disorder is induced 
by the input laser beam), the behavior of light can be in- 
terpreted using the same paradigms of modern statistical 
physics. 

If a light filament, or spatial optical soliton, can be 
treated as a classical particle, many interacting spatial 
solitons correspond to a liquid or high density gas. If 
some noise is present in the system, a Brownian dynam- 
ics model is readily introduced, and the filaments behave 
like particles dispersed into a solvent, which is one of the 
simplest definitions of soft-matter. The dynamic phase 
transition is the natural way to describe phenomenolog- 
ical transformations, like the formation of clusters. This 
has been shown by numerical experiments in this article 
and experimentally in future publications. 

Nonlinear optics can be hence used to test, theoreti- 
cally and experimentally, some of the ideas of the physics 
of complexity, as well as the latter can be used to explain 
many high-field phenomena, like laser-filaments genera- 
tion and interaction. In this manuscript these ideas have 
been applied to explain a possible manifestation of com- 
plex light. 
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